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ABSTRACT 

The Hydra iV-body hydrodynamics code has been modified to model, from 
the end of the Sedov phase, the effects of supernovae on the surrounding 
medium. The motivation is to investigate the feedback of energy into the 
interstellar/intergalactic medium. We compare our results for supernova 
remnants (SNRs) in a uniform medium to previous detailed work on the late 
evolution of SNRs. The code is found to reproduce the bulk characteristics of 
SNRs well. Results on the effects of a single central SNR on Plummer clouds are 
presented. The feedback of kinetic energy and the percentage mass loss can be 
parameterised in terms of the cloud mass and characteristic radius in a simple 
way. 

The kinetic energy fraction returned to the ISM from a SNR is < 3 per cent. 
The removal of gas from the cold, dense phase and the addition of energy due 
to the lowering of the potential energy of a cloud is at least as significant, if not 
much more so, than the kinetic energy leaving a cloud. 

Subject headings: supernova remnants — supernovae:general — hydrodynamics 
- methods: numerical 
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1. INTRODUCTION 

The mechanism by which energy and metals are fed back into the inter-stellar medium 
(ISM) by the supernovae of massive stars is an interesting and complex problem with 
applications to many areas of astrophysics. In this paper we present an adaptation of the 
Hydra iV-body hydrodynamics code (Couchman, Thomas & Pearce 1995) which is capable 
of modelling the effects of multiple interacting supernovae over a large range of scales. 

Hydrodynamical simulations of galaxy evolution and cosmological structure formation 
have recently been a topic of particular interest (eg. Katz 1992; Katz, Hernquist & Weinberg 
1992, 1999; Navarro & White 1993; Steinmetz & Muller 1995; Frenk et al. 1996; Pearce et 
al. 1999). Feedback from supernovae appears to be required to reheat gas and prevent the 
'cooling catastrophe' (White & Rees 1978; White & Frenk 1991), caused because cooling 
is very efficient in small halos where the cooling time is less than the dynamical time. As 
the particle mass in a simulation increases, more and more of these halos are resolved and 
a greater fraction of the available gas cools, in contradiction with the observation that 
only a small fraction of baryonic material is in the form of stars or cold gas. Presently the 
efficiency at which energy is returned to the ISM and beyond to the intergalactic medium 
(IGM) is a free parameter with little more than handwaving arguments to support the 
values adopted (eg. Katz 1992; Navarro & White 1993; Mihos & Hernquist 1994; Gerritsen 
& Icke 1997). It is to be expected that such efficiency parameters depend upon the density 
and metallicity of the surrounding gas at the very least. One of the main aims of this 
project is to better constrain the values of these feedback parameters (the amount of energy 
and material returned to the ISM as well as the spatial extent and form of the feedback) 
for a wide range of stellar systems. 

An understanding of the early evolution of star clusters requires knowledge of how 
the residual gas, left within a cluster after star formation has occurred, is expelled (eg. 
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Tenorio-Tagle et al. 1986; Goodwin 1997a, b). Such an understanding would also help 
provide constraints on the initial conditions and formation models of star clusters (Goodwin 
1997b). Within a larger context, information on the effects of multiple supernovae is 
required to construct detailed stimulated star formation models. 

Much detailed work (both hydrodynamic simulations and analytic calculations) has 
been carried out upon the dynamics of individual supernovae expanding into a uniform ISM 
(see Lozinskaya 1992 and references therein), however, little of this work has concentrated 
upon the late phases of that expansion (exceptions including Cioffi, McKee & Bertschinger 
1988; Slavin & Cox 1992; Thornton et al. 1998). Even less has been published on the 
effects of supernovae in gas clouds; analytic calculations were made of the effects of multiple 
supernovae on gas clouds by Dopita & Smith (1986) and Morgan & Lake (1989), but due to 
the complexity of the problem the treatment was understandably simplistic. More recently 
Petruk (1999a, b) has published simulations of the very early evolution of a SNR in a 
density gradient. 

In this paper we describe the extensions we have made to the Hydra code to model 
this problem and present a number of convergence tests that we have performed. We also 
compare the results from this code to previous analytic and simulation work on individual 
supernovae. In section 2 we describe the testing of the code and show that it replicates the 
behaviour of a supernova remnant expanding into a uniform medium as studied by other 
authors. In section 3 we study the effects of a central supernova in Plummer clouds of 
various masses and characteristic radii. In section 4 conclusions are drawn and the future 
applications of the code detailed. 
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2. TESTING 

In this section we present the code we have used and test both its self-consistency and 
it's ability to reproduce known results for the evolution of SNRs in a uniform density ISM. 
Firstly we outline the modifications to the Hydra code. Then we show that the evolution 
of the energy content of an SNR and its expansion law is in good agreement with other 
simulations. We show that these laws hold for different ambient densities and that the 
properties of an SNR scale with its initial radius. 

We do this to show that we are confident of the code's ability to model the evolution 
of a SNR when applied to the problems of SNR evolution in a cloud in Section 3 (and in 
future papers). 

2.1. Hydra 

The public release version of Hydra, upon which the code used in this work is based, 
has been extensively tested in a variety of astrophysical fields (Couchman 1991, Couchman, 
Thomas & Pearce 1995, Thacker et al. 1999, Kay et al. 1999, Benson et al. 1999). 
It is an adaptive particle-particle, particle-mesh code that includes smoothed particle 
hydrodynamics (SPH) (Gingold & Monaghan 1977; Lucy 1977) to follow the gaseous 
component. 

For this work we have chosen to use the latest public release version of Hydra, 
(Hydra3.0). This version incorporates a standard pairwise Monaghan type r.v viscosity 
rather than one based on the divergence of the local flow as used in earlier release versions 
of Hydra, because the r.v viscosity provides better shock capturing. We also employ the 
improved neighbour counting estimator as described by Thacker et al. (1999). To follow 
the radiative cooling of the gas we use the cooling rates calculated by Sutherland & Dopita 
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(1993) interpolated to the chosen gas metallicity. 

In addition we have rewritten the SPH algorithm in order to perform complete 
neighbour counting even in low density regions. Previously, neighbour searching within 
Hydra was limited by the search length imposed by the size of the particle-particle grid, 
resulting in an uncomfortably high minimum resolved density. For our problem we wish to 
resolve the hot, low density gas within expanding supernovae shells and so have removed 
the maximum search-length restriction. 



2.2. Individual supernovae 

Here the mechanism for adding individual SNRs into a simulation is outlined. 

The initial stages of a SNRs evolution are complex as an enriched, very high 
temperature gas ploughs into the surrounding ISM. We do not wish to simulate this phase 
of SNR evolution. In order to deal with multiple SNRs over a large range of length scales 
we ignore evolution prior to the onset of the late stages of evolution characterised by the 
pressure-driven snowplough (PDS) phase (cf Cioffi et al. 1988). 

The choice of the beginning of the PDS phase as the starting point of the evolution 
within the code is due to a number of considerations. Firstly we wish to simulate the 
evolution of many SNRs in complex environments - the mass and time resolution (as well as 
the complex input physics) required to simulate the earlier stages of evolution are beyond 
our remit. Secondly, the evolution of a SNR to the PDS phase is reasonably well understood 
and provides a good basis for our simulations, allowing us reduce the mass resolution 
required to model its evolution. Lastly, the evolution of a SNR through the adiabatic Sedov 
phase even in a strong density contrast such as would be found in a molecular cloud is 
expected to be roughly spherical; significant non-spherical evolution would only be expected 
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in the late PDS phase (Dopita & Smith 1986). In dense environments, such as molecular 
clouds, the radius at which the PDS phase starts is much less than 1 pc, thus for significant 
non-spherical evolution during the Sedov phase the density change must be on a similarly 
small scale. 

The early evolution of SNRs comprise a short-lived initial free-expansion phase followed 
by an adiabatic expansion (the Sedov phase) once the mass of swept-up material exceeds 
the initial ejecta mass (Spitzer 1978). We will concern ourselves only with the later stages 
of SNR evolution after the Sedov phase has finished. 

A few xlO 4 yrs after the initial supernova, the temperature of the adiabatically 
expanding remnant falls to the point at which cooling is efficient in the outer regions: 
T (5 — 6) x 10 5 K. At this point around half of the remnant mass forms a thin shell which 
expands into the ISM in a pressure-driven snowplough (McKee & Ostriker 1977; Lozinskaya 
1992) . It is at the point at which the PDS phase begins that we will start to simulate the 
evolution of SNRs. 

At the beginning of the PDS phase the thermal energy E? in the remnant is 
E T PS 0.36E 51 . The temperature T PD g of the interior gas is given by 



T PDS = 1.545 x 10 10 ^i^ K, (1) 
«o*p DS 

typically a few xlO 6 K. This leads to a value of the interior pressure in very good agreement 
with that of Chevalier (1974). The thin shell (containing half the mass) is at around 5 x 10 5 
K. 

Initially a supernova inputs energy into the ISM such that the thermal energy is 
Et ~ 0.72_E SN and the remainder is in the form of kinetic energy where E S n is the total 
supernova energy, where the usual value is -Esn = 10 51 ergs (Chevalier 1974). By the time 
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a remnant enters the PDS phase, around 50 per cent of the thermal energy will have been 
lost (Lozinskaya 1992). The shell temperature is (5 — 6) x 10 5 K and the interior pressure 
will have fallen to P ?a .E SN /(47ri?p DS ) where rpos is the radius at which the transfer to the 
PDS phase occurs (McKee & Ostriker 1977). 

The radius r PDS at which the PDS phase starts is given by Cioffi et al. (1988) as 

_ 1 A Z?0.29^, -0.43 A-0.143^_ (n\ 

r PDS - 14£ SN n Q m pc (2) 

where Uq is the ambient density in Hydrogen atoms cm -3 and ( m is related to the metallicity 
(for Solar metallicity ( m = 1). 

The initial velocity vpds of the shell is given by 

v PDS = 413^ 43 £ 5 7 n Cr 14 km s" 1 (3) 

These parameters (radius and velocity) are similar to those given by other authors (eg. 
Chevalier 1974; Falle 1981; Blondin et al. 1998). Differences are mainly due to the use of 
different cooling functions although the shell velocity predicted at the same radius is very 
similar for all calculations. 

The main reasons for preferring the calculations of Cioffi et al. (1988) over other 
calculations are that they include a correction factor for non-Solar metallicities and their 
good agreement with the simulations of Thornton et al. (1998). 

2.3. SNRs in the code 

In the remainder of this section we compare our results for a SNR in a uniform medium 
with those of other authors as well as testing the convergence and self-consistency of our 
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results. 

SNRs are set-up by taking all of the particles in a sphere of radius tpds centred on the 
position of the supernova. Half of the particles are unmoved but are heated up to Tpds to 
represent the hot interior (uniform density and pressure, no bulk motions). The other half 
of the particles are distributed on the surface of a sphere of radius rpos and given a velocity 
outwards from the centre of the sphere of upos at a temperature of 5 x 10 5 K. 

Low particle numbers within rpos have an effect on the evolution of the SNR due 
to shot noise causing significant non-sphericity. Generally more than « 60 particles are 
required within rp^s to model the SNR well although as few as 20 do a reasonable job (we 
would not wish however to trust simulations where the included particle number is this 
low). 

2.3.1. SNR energy 

Table [3] gives the parameters used in 10 example runs covering a factor of 8 in N, 2.5 
in box size and over 30 in mass resolution for the same physical initial conditions. These 
runs are used in the next 2 subsections to illustrate the stability of the code at reproducing 
results over this range of parameters and to compare with other authors to show the ability 
of the code to well represent the post-Sedov phase evolution of a SNR. 

Figure [l] shows the change in the total energy within the simulation volume, normalised 
to the value at 1 Myr to take into account the different amounts of total thermal energy 
present in simulation volumes of varying size (in large boxes the thermal energy of the 
undisturbed gas is by far the largest contributor to the total energy). The 'lost' energy 
has been radiated away by the hot gas, the energy conservation being good, a maximum of 
0.1 — 0.2 per cent. Good agreement is seen over a wide variety of particle mass and box size. 
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The small differences are mainly due to the requirement that there be an integer number 
of particles in the interior and shell of the SNR. This leads to small differences between 
the initial thermal and kinetic energies of each model which amount to a few per cent, a 
difference which is amplified as time progresses. For instance, the initial kinetic energy 
of runs 1019 and 1020 differ by 4 per cent, by 6 Myr the difference in fig. [I] of 0.02£ , 51 
represents approximately 4 per cent of the total energy change. 

A more detailed examination of the energy balance for three of the SNRs is shown 
in fig. H, approximately 0.6 Myr after the supernova explosion. The solid lines show the 
sum of the thermal and kinetic energy within each radius with the dashed and dot-dashed 
lines showing the total thermal and kinetic energy respectively within each radius. The 
outwardly moving shell is clearly visible. The differences in interior thermal energy are due 
to shot noise and very low particle numbers in these regions - the high thermal energies of 
runs 1002 and 1017 within the inner 10 pc are due to 1 hot particle. The agreement of the 
energies over a factor of more than 30 in mass resolution is good and is representative of all 
runs. The shell has a peak density at approximately 60 pc at 0.6 Myr and is wider where 
the mass resolution is poor. 



2.3.2. SNR expansion law 

The expansion of a SNR can be roughly characterised as a power law of the form r oc t v 
(with t being measured from the time of the supernova). McKee & Ostriker (1977) find for 
negligible external pressure that v = 2/7 (~ 0.286). Chevalier (1974) finds v = 0.305 for 
early times and for late evolution (t > 0.75 Myr) v = 0.32. This is close to the Cioffi et al. 
(1988) offset solution where r oc {t — t ffset) ' 3 - From inspection of fig. 3 in Thornton et al. 
(1998), they appear to obtain a value around v = 0.3. All these values are below the power 
law for a purely momentum conserving snowplough where v = 0.4 (Spitzer 1978). 



- 11 - 



Figure |3| shows the evolution of the shell radius with time for the ten convergence test 
runs (which model the same physical conditions at different mass resolutions and box sizes) 
detailed in table [l| In fig. [| the shell radius is determined by the mean radial distance from 
the centre of the SNR of the densest particles. At low shell radii r this is less reliable due to 
the smaller number of particles in the shock (especially when the mass resolution is poor). 

The best fit to the first 1.5 Myr of evolution is shown by the dashed curve fitted with; 

r(t)=r PDS (^- + l)" (4) 

with r in pc and t in Myr (t is the simulation time, ie. the time since the onset of the PDS 
phase), tpDs is of order the duration of the PDS phase. Leaving tpds, ^pds and v as free 
parameters the best fit is given by rp^s = 18.34 pc, tpDs = 0.012 Myr and v = 0.317. For 
an n = 0.5 cm -3 ISM, r PDS is actually 18.85 pc (from eqn. and ^ds = ^pds/^pds is 
around 0.02 Myr (Ciofli et al. 1988). 

The late stages (> 1.5 Myr) are well fitted by the solid line, r = 55.8 + 16. 8t, a 
reasonably good fit to free expansion of a sound wave, the sound speed in this case of 15.6 
pc Myr -1 . 



2.3.3. Different ambient densities 

We have investigated the differences in SNR evolution within ISMs of different ambient 
densities. Table ^| shows the range of ambient densities tested, covering 4 orders of 
magnitude. Previous studies (with the exception of Thornton et al. 1998) have concentrated 
on investigating low ambient densities representative of the ISM in the Solar neighbourhood. 

Figure |4] shows the shell evolution for the range of densities presented in table 0. Also 
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shown is the radius and time at which the shell velocity falls to the sound speed for each 
density. 

From the data presented in table ^ the radius at which the sound speed is reached is 
fitted by a power-law in density of r sound = 65uq °' 37 pc and the time by a power-law of the 
form t sound = 1.2no a38 Myr. 

These power-law fits to r soun d and t SOU nd are almost exactly what are expected. 
Using eqn. £| and the known dependencies on n of rpDs(oc tiq 0A3 - eqn. §) and 
^pds(oc rpos/wpDS oc n - 57 -eqns. | & D gives for t/t PDS » 1 

r oc n - a25 t a32 (5) 

as v soun d = r soun d/t SO und = constant leads to the predictions that t SOU nd oc n^ 036 and 
^sound oc ng °' 37 , very close to the fitted relations from the simulations. 

The initial energy (thermal and kinetic) of the SNR and the evolution of the total 
energy is approximately the same for each density if the evolution of a SNR through the 
PDS phase is scaled so that the radial properties are in units of rp^s as illustrated in 
fig. |5|. The radial total, thermal and kinetic energies are plotted against r/rpDs for 3 models 
spanning a range of 4 orders of magnitude in density for the time when r s h e ii ~2x tpds- 
For n = 0.01 cm -3 , 2rpDs ~ 200 pc, for n = 0.5 cm~ 3 , 2rpDs ~ 38 pc while for n = 100 
cm -3 , 2r PDS w 3.9 pc. 

As is clear in fig. [5] the energy profiles of the SNRs in terms of r/rppjg are very similar, 
indeed the main differences are due to the fact that the output radii shown are not exactly 
r s heii = 2 x rpDs but differ by a few per cent (due to the finite number of outputs). 

When r = 2 x tpds the mass of the hot interior is around twice the interior mass at 
the onset of the PDS phase as most of the swept-up material remains in the dense shell. 
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The interior temperature falls by a factor of around 2.7 while the shell temperature has 
fallen to ~ 10 4 K. The shell velocity has fallen by a factor of 5 (as expected if v scales as 
^pds^Apds + 1)~ ' 68 , cf. eqn. f|). The scaling of the SNR properties with tpds across all 
densities is because the energy content of SNRs is density independent - a doubling of the 
volume of the SNR will result in the same energy changes throughout the SNR independent 
of the ambient density (as swept-up mass dominates). 

2.3.4. Late-starting SNRs 

Figure |^ shows the late evolution of three SNRs, two started at the onset of the PDS 
phase and the other started at r = 2 x tpds with T s h e ii = 10 4 K and the interior temperature 
given by T SN r = 0.37 x T PDS (where T PDS is given by eqn. |) and v sheU = 0.2v PDS . The 
total energy within the late-starting SNR is very similar to that in those that started at 
rpDs but the relative contributions of the kinetic and thermal energies within the bubble 
are slightly different. This is due to the shorter amount of time which the 2rpDs simulation 
has had to convert bulk motions into thermal energy within the bubble. Again the main 
differences between the simulations are due to outputs not being at exactly the right times. 
The ability to start modelling each SNR at 2 x rp^s will become very useful when we wish 
to add SNRs to massive gas clouds reducing the required mass resolutions in simulations by 
a factor of 8 (see section 3). 

2.3.5. Metallicity effects 

The effect of different metallicity environments can be significant during the early 
evolution of SNRs, during the PDS phase however the differences due to metallicity are 
found to be negligible (see also Thornton et al. 1998). Initially a metallicity of ( m = 0.01 
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will make a factor of 1.93 difference in r PDS and 0.37 in w PDS as compared to a ( m = 1 
SNR (from eqns. 2 & 5 respectively). However, the effect of these different initial SNR 
conditions is not as significant as it might first appear, because the velocity of a ( m = 1 
SNR at r = 1.93 x r PD s(Cm = 1) is « 0.22t> PD s(Cm = 1) (see above) resulting in only slightly 
different expansions. 

The subsequent evolution of SNRs with different metallicities through the PDS phase 
are virtually indistinguishable (as also found by Thornton et al. 1998). The evolution 
of the shell radii is slightly different at early times due to different values of r PPJ s but as 
mentioned above this difference is soon nearly cancelled out and at late times the shell 
radii differ by only a few pc. When SNRs are started with the same r PD s but different 
metallicities the difference is negligible. The evolution of the total energies within each SNR 
are also remarkably similar. Differences in metallicity affect the cooling of the hot interior 
and hence the pressure and the extra driving force on the shell. However, during expansion 
into a uniform ambient medium the shock acceleration is dominated by the deceleration 
due to the sweeping-up of material, making metallicity (which effects the pressure-driven 
acceleration) a minor factor in the late-time evolution of the SNR. 

The differences in energy between different metallicities are accounted for by the 
different initial conditions at the onset of the PDS phase however it is not clear that even 
the early stages of SNR evolution would be affected by low-metallicity effects. A type II 
supernova with a progenitor mass of 25M Q will produce around 2.4M of heavy metals 
(Tsujimoto et al. 1995), the majority (1.8M ) being O 16 . In the case of a SNR expanding 
into an ISM of density n = 100 cm -3 (where the mass within r PD s is ~ 75.5M ) this 
will enrich the swept-up material from ( m = 1 (Solar) to ( m = 2.6 and from ( m = 0.01 to 
Cm = 1.6. Such enrichment is very significant and will alter the evolution of the pre-PDS 
phases. However, by 5r PDS the amount of swept-up material is 944OM and has been 
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enriched from ( m = 0.01 to ( m = 0.023, a fairly insignificant (from the point of view of 
cooling) amount. Thornton et al. (1998) did not find a significant effect when including 
the metals ejected from the supernova (Thornton 1999). Obviously this enrichment and its 
history is of vital importance for understanding chemical evolution, and the enrichment 
history can be followed by this code (once a number of assumptions about mixing have 
been included), but it has little bearing on the late-time evolution of an individual SNR in 
a uniform medium. 

2.3.6. Computational aspects 

The softening included in the simulations does not greatly affect the results of the 
shell expansion and energy transfer. Similar runs differing in softening are virtually 
indistinguishable in their results (two such runs, differing by a factor of over 4 in softening 
- 1002 & 1005 - are included in fig. [l]). The softening selected for simulations should be 
as small as is reasonably practicable given that the smallest region in which results can be 
believed in detail is no less than the softening length, which also determines the width of 
any shock fronts that are present. Unfortunately the number of timesteps required rises as 
the softening length is decreased. 

The width of the shell that forms the snowplough is difficult to determine accurately as 
it becomes slightly aspherical as it evolves due to initial Poisson fluctuations in the particle 
positions. This effect is worse for the models which began with only a few particles in the 
SNR. When spherically averaged profiles are produced this leads to an artificial broadening 
of the snowplough and general smoothing of the steep shock fronts. The snowplough width 
is expected to be a few pc (Thornton et al. 1998), in these simulations the shock width is 
approximately the softening. 
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2.4. Summary 

In this section we have shown that the code is stable and converges for reasonable 
selections of box size and particle mass. The code is able to reproduce the results of other 
authors for the situation of a SNR expanding into a uniform ambient medium during the 
PDS phase of evolution. In the next section we will place SNRs in the centres of Plummer 
model gas clouds. 



3. CENTRAL SUPERNOVAE IN A GAS CLOUD 

While the evolution of a SNR in a uniform density ISM provides a good test of the 
ability of the code to reproduce the characteristics of a SNR as found in more detailed 
hydrodynamic calculations, it is not physically realistic. The effects of a supernovae in a 
stratified medium (eg. a gas cloud) is a more practical problem. Massive stars are found in 
young clusters, their lifetimes being so short that they are expected to still be embedded in 
gas remaining from star formation. The significant evolution of a young SNR is therefore 
expected to occur in the confines of its parent GMC. Massive stars are also found to be 
centrally-concentrated in clusters (Hillenbrand 1997; Carpenter et al. 1997) so that an 
investigation of a central SNR is a good approximation to the evolution of the first SNR in 
a cluster. 

Some analytic calculations of the effects of supernovae in a gas cloud have been made 
by Dopita & Smith (1986) and Morgan & Lake (1989) to estimate the number of supernovae 
required to totally disrupt a gas cloud of a given mass. Morgan & Lake (1989) using more 
detailed cooling functions than Dopita & Smith (1986) found that the minimum mass of a 
(1/r 2 ) cloud for it to confine a single central supernova was 4 x 10 4 M Q . 

As noted in section 2.3, Dopita & Smith (1986) find that a SNR within the high-density 
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environment of a gas cloud will remain roughly spherical during the adiabatic phase of its 
evolution. Using this result we may place SNRs in the cloud at the start of the PDS phase 
in the same way as they are placed in the uniform medium simulations. The high densities 
in gas clouds (no > 10 4 cm -3 ) means that the adiabatic phase will end whilst the shell 
radius is very small (<C pc). 

Cloud initial conditions are based on a Plummer model (n = 5 polytrope) with mass 
distribution; 

Mr 3 1 

where M is the total mass and R is the scale length. The half-mass radius of a cloud is 
then « 1.3-R. Clouds are constructed out to a maximum radius r max of a few R (normally 
r max = 20 pc), resulting in the actual mass being slightly less than M. 

In this paper we will only deal with Plummer model clouds. Obviously the effects of a 
SNR will depend, at least to some extent and possibly very significantly, upon the density 
distribution of the parent cloud. These will be dealt with in detail in a paper to follow. 

Observationally, Giant Molecular Clouds (GMCs) can approximated very roughly by 
clouds with M = a few xlO 4 to a few xlO 5 M and R — a few pc (eg. Harris & Pudritz 
1994) which we used as initial conditions for our simulations (full details in table 0). Before 
inserting any supernovae, our clouds are allowed to relax to a stable state in which they are 
virialised. Support is provided by the thermal energy of the gas and the bulk kinetic energy 
is negligible (although the temperature of the gas represents the turbulent velocity which 
supports the cloud rather than the molecular temperature). The clouds are self-gravitating 
as are real giant molecular clouds (Rivolo & Solomon 1988) and so require no pressure 
confinement from an external hot diffuse gas. 
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The insertion of a supernova has the effect of adding considerable amounts of kinetic 
and thermal energy to the cloud. By the end of the Sedov phase, the extra energy amounts 
to around 6.4 x 10 50 ergs, comparable to or greater than the potential energy of the system 
(roughly 8.6 x lO 4O (M/M ) 2 /(i?/pc) ergs). The shock from the supernova then passes 
through the the cloud. 



3.1. Some numerical considerations 

The number of particles within tpds required for the SNR evolution to be a good 
approximation to more complex simulations is > 40 with > 60 giving the best results (see 
section 2). For a Plummer model cloud the central density is 



3M 

Po = 4^ (7) 

where M is the mass of the cloud and R is the Plummer scale length. Using equation ||] the 
mass interior to r PDS at the cloud core will be 



R 6 / 7 

M PDS « 5.48^ M (8) 
M 5 2/7 

where R pc = R/pc and M 5 = M/1O 5 M (we shall be using R pc and M5 throughout the rest 
of this section). We require a mass resolution at least 60 times smaller than this giving, for 
R pc = 3.5 and M 5 = 1 a minimum particle number of N > 3.74 x 10 5 . 

Such resolutions, while attainable, require significant computing time to run. Our aim 
here is to introduce SNRs in such a way as to allow many simulations covering a variety 
of initial conditions to be completed, so as to explore the parameter space. Starting SNR 
evolution at r = 2 x tpds (see Section 2.3.4) allows a factor of 8 reduction in the mass 
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resolution and therefore the number of particles required to model the cloud. In the case of 
a 10 5 M Q cloud 2 x 32 3 particles are required, on a workstation each timestep takes around 
3 minutes resulting in a 20 Myr simulation in around 30 hours. 

3.2. Feedback as a function of cloud structure 

Here we present the analysis of the energy and mass feedback as a function of M 5 
and R pc , the cloud models are presented in table |3|. Initially all clouds are virialised and 
particles are said to be lost once they cross the box edges which are a distance 2r max from 
the centre of the cloud (at a distance roughly equal to the tidal radius of a GMC). It is 
these lost particles that comprise the feedback from the cloud into the larger ISM. In most 
cases the SNRs are started at 2rpDs- 

As a 'standard' model we use the run with parameters M 5 = 1, _R pc = 3.5 to 
illustrate the general features of SNR evolution in a Plummer cloud. Initially the cloud has 
0.0356 x 10 51 ergs of thermal energy and a potential energy of —0.0711 x 10 51 ergs. The 
SNR starting at 2r PDS adds a total of 0.123 x 10 51 ergs of kinetic energy and 0.289 x 10 51 
ergs of thermal energy (calculated from section 2) creating a net positive energy for the 
system of 0.38 x 10 51 ergs. 

3.2.1. The evolution of the SNR 

The evolution of the shell velocity with radius for the standard model is shown in 
fig. 0. It is typical of all runs that the velocity drops from its initial value (wpos or 0.2t>PDs 
of started at 2rpDs) reaching a minimum value f m ; n at around 2 to 3 R pc before accelerating 
to around 1.5 x t> min on leaving the box. 
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It is found later that t> m in and associated quantities may be used to parameterise all 
of the effects of a SNR on a gas cloud. In this subsection we sketch a model of the SNR 
evolution that predicts t> m in in terms of the two cloud parameters M 5 and R pc . 

The velocity reaches v m { n when the acceleration due to the interior pressure matches 
the deceleration due to the sweeping-up of material. The force on the shell due to pressure 
F press (r) goes as 



F press = (47rr 2 )P (9) 

where P is the pressure when the shell radius is r(t). 

The initial radius of the shell (Vpog) is very small compared to R and as adiabatic 
cooling goes as (rpDs/^) 5 the majority of the cooling will be adiabatic, especially at late 
times so 



P - ^PDS^f (10) 



with P given by p kT PDS (where p is the central density oc M 5 /R^ c ). The initial 
temperature Tpds is given by eqn. [I] as 



M 0.29 

T PBSOC^ (11) 
-"'pc 

The force due to the sweeping-up of material (momentum conservation) F mom (r) is 



(r) = (Airr 2 ) p(r)v (r) 



(12) 



equating eqns. |9| and [12] using eqn. [10| gives 
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P(^in)< m = Po4^ (13) 



/v>5 
min 



using eqn. |Tl| and assuming that r min 3> R gives 



^min ^ ^0.93 (^J 

Fitting the powers of M 5 and i? pc simultaneously gives a best linear fit to t> m in of 



R 032 

6-58^ km s- 1 (15) 



which is compared to the simulations in fig. 0. The similarities between eqns. 14 and 15 



are remarkable considering the very simple assumptions that went into the formulation of 



eqn. 14 



The accelerating force due to the pressure goes as r -5 and so rapidly becomes negligible 
at high r. It is still able, however, to increase the velocity of the shell by the time it leaves 
the cloud to ~ 1.5i> mm . 

The above rather simplistic treatment avoids (as is necessary) many details of the 
actual SNR evolution in a cloud. In the simulations the main shell of the SNR is followed 
by a weaker shock which is created by the shocking of infalling material filling the interior 
when it reaches the centre (and meets other infalling material). This creates a very complex 
density /pressure structure in the cloud as the SNR evolves. 

We know go on to show how the effects of a SNR on a cloud can be parameterised very 
simply in terms of iVn and closely related quantities. 
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3.2.2. Metallicity effects 



As the temperature factor T is so important in setting the value of t> m in, then 
metallicity, Z, is an important factor in that in clouds of lower Z than the fiducial Solar 
Z — 1, cooling will be less efficient and T and t> m m higher. In the case of a Z = 0.01 
metallicity cloud the kinetic energy lost increases to 2.5 x 10 49 ergs and the mass lost to 
5.4 x 1O 4 M , a 67 per cent increase on the Z = 1 cloud in both cases. For a Z = 0.1 cloud 
the increase in both values is 46 per cent. This would indicate that the effects of the first 
generation of supernovae would have been much more dramatic than those which occurred 
later once metal enrichment had taken place. 



The mass loss from the cloud would be expected to be related the escape velocity from 
the cloud at the point where the minimum velocity is reached. The escape velocity at 2R pc 



is given by 18. 6^/ M 5 /R pc km s -1 and it is found that if t> m in > f C sc(?"mm) then the cloud 

completely disrupts, ie. all of the mass is lost. At the other extreme if v m i n < 0.2v esc (r min ) 
then no mass loss (and hence no feedback) occurs. In the intermediate regime the mass loss 
is related to v esc (r min ) by 



3.2.3. Mass loss and disruption 





(16) 



as illustrated in fig. |T(]. 
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3.2.4- Energy feedback 

The amount of energy returned to the ISM from a cloud would be expected to be 
related to the kinetic energy of the cloud at t> min given by T min = l/2M she nf where M shell 
is the mass of the shell at r min , most simply we might expect M shell oc M 5 . 

Illustrated in fig. |TT] is the feedback energy from table [| against T min /10 51 ergs 
= 0.001M5(f m i n /km s -1 ) 2 . Two linear regimes are obvious in the behaviour of the feedback 
energy In the upper limit where the cloud has totally disrupted the feedback energy is 
fitted by 



E lost = -0.0012 + 0.25M 5 ^ in x 10 51 ergs (17) 
in the lower limit where the mass loss is negligible the feedback energy is 

E elost = -0.012 + 0.25M 5 < in x 10 51 ergs (18) 

as can be easily seen the slope of these two relationships is the same and the constant 
factor, different by exactly one order of magnitude, is the only difference. The intermediate 
stage marks the change from the total cloud destruction and high feedback to the low mass 
loss and low feedback regime. This region is only a very small region of the total M5 — R pc 
parameter space. 



3.2.5. The final state of clouds 

As stated above when t> min > t> esc a single SNR is enough to disrupt a cloud completely 
and when v m \ n < 0.2t> esc no mass loss occurs at all. In the intermediate regime some gas is 
retained in a bound object - a new cloud with lower mass than the original. 
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The ratio of the final to initial potential energy, flf/fli is related again to the ratio of 
the minimum velocity of the shell and the escape velocity v min /v csc 

^ = 1.26 - 1.05— (19) 

* 'i ^esc 

as illustrated in fig. This relationship is only valid for 0.3 < v min /v csc < 1, beyond 1 the 
cloud is totally destroyed and Qf/Qi = and less than 0.3 flf/fli asymptotes to 0. 

The change in potential energy of a gas cloud due to a central supernova can sometimes 
be significantly more important than the feedback of kinetic energy as illustrated in table [| 
Clouds lose all of their initial potential energy if they are completely disrupted. If the cloud 
loses only part of its mass then the significant energy change occurs as a loss of potential 
energy from the cloud (a net gain of energy). Even when no mass loss or feedback occurs a 
small change (of the order of a few percent) in the potential energy is observed as the cloud 
expands slightly. 



As fig. [13] shows, the remaining material settles back into a new configuration that can 
also be described as a Plummer model with a lower central density and larger scale length. 
Thus knowing the initial mass M and characteristic radius R gives the final characteristic 
radius Rf via; 

Knowing the final mass and characteristic radius of the cloud allows the effects of further 
central supernovae on that cloud to be calculated - as long as the interval between 
supernovae is greater than the time required for the cloud to relax to a new equilibrium. 

The time taken for a 5 x 1O 5 M cloud to recover from a central supernova is very short, 
the majority of the SNR's kinetic and thermal energy is radiated away in < lMyr. Assuming 
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a Salpeter IMF, the number of massive stars within a cloud that will go supernovae is 
-^SNe ~ 0.006Me where M is the mass of the cloud and e is the star formation efficiency. 
Taking e = 0.01 and M = 5 x 1O 5 M gives iVsNe = 30. If these supernovae occur centrally 
and evenly spaced over a 30 Myr period then no ejecta will escape the cloud and the cloud 
will not be disrupted by the SNRs. 

These results hold for Plummer model clouds that can be characterised solely by M 
and R. This study will be extended to other density distributions in a future paper which 
we are preparing. 

4. CONCLUSIONS 

The Hydra iV-body SPH code has been extended to allow the simulation of the 
evolution of supernova remnants (SNRs) from the onset of the pressure-driven snowplough 
(PDS) phase. In section 2 this code was seen to be able to produce convergent results over 
a wide range of parameter space and reproduce the results on the evolution of SNRs from a 
variety of previous authors. The power of this code is the ability to simulate the evolution 
of SNRs in a variety of environments using a workstation in a reasonable time (of the order 
of days). 

This code represents the first time that we are able to model multiple, interacting SNRs 
in gas clouds and complexes. This paper is the first detailed analysis of SNR evolution 
in a gas cloud to examine the feedback parameters so important in galaxy formation and 
evolution calculations. 

In section 3 we presented new results on the effect of a single, central supernova 
on Plummer model gas clouds of various masses and characteristic radii. The results of 
Dopita & Smith (1986) and Morgan & Lake (1989) were found to be too simplistic. The 
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evolution of a central SNR is very complex, the late stages of evolution are governed by 
pressure-driving from a hot interior accelerating the SNR shell down the density gradient. 

Feedback from a cloud is defined to be the mass and energy which pass out of our 
simulation box (whose size is approximately the tidal radius of the cloud). The results may 
be summarised as: 

• The efficiency of energy feedback, mass loss and cloud destruction for a central 
supernova in a Plummer cloud of mass M and characteristic radius R is related to the 
minimum velocity i> m in that the shock reaches during it's evolution. This minimum 
velocity is given by v mhl = R^f/M^ 77 . 

• For t> min > t> esc the cloud is totally destroyed, while for t> min < 0.2t> esc the SNR is 
completely contained and no feedback occurs. 

• The mass lost from a cloud is related to the ratio of the minimum velocity to the 
escape velocity twAw as M Xost /M = 0.25(v mhl /v csc ) 2 - 5 . 

• The energy feedback £i ost has two main regimes with equal slopes where 

E lost = -C E + 0.25M 5 V^ in x 10 51 ergs where C E = 0.0012 when v min /v csc > 1.1 and 
Ce = 0.012 when f m i n /f e sc < 0.9. 

• The loss of (negative) potential energy is often the largest (positive) increase of energy 
in the system and at least of order the feedback of kinetic energy, the feedback of 
thermal energy being negligible in comparison to both. 

• The final state of a cloud that is not destroyed is close to a Plummer 
model with final characteristic radius Rf related to the initial parameters by 

R f = (i - M lost /M) 2 R(n f /n t )-\ 
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• The efficiency of feedback increases rapidly with decreasing metallicity suggesting 
that feedback at early epochs was far more efficient. 

Our simulations are obviously not perfect. The code takes no account of magnetic 
fields (following standard astrophysical practice) which may well be important, especially at 
late times. In addition the feedback from massive stars into the cloud before they become 
supernovae is neglected, although in the dense environments of cloud cores we may be 
justified in ignoring this effect (Franco, Garcia-Segura & Plewa 1996). A simple method of 
including such effects in the code is being developed. Despite this we believe that this code 
represents a significant step in modelling the effects of multiple SNRs on clouds and the 
larger ISM. 

4.1. Future directions 

As noted in the introduction, a code such as this has many interesting applications. 
It can be used to investigate how gas clouds of different sizes and shapes are affected by 
internal supernovae, calculating the feedback of energy and mass into the ICM. We are in 
the process of preparing papers that will expand the current work into an investigation of 
off-centre SNRs in gas clouds and of the effect of different density distributions on feedback 
parameters. 

In addition we will be able to investigate how the cloud is disrupted and in what way 
gas is expelled on a small scale. This has important consequences for investigating the early 
evolution of star clusters. A paper is in preparation on the effects of gas expulsion on the 
dynamics of the stellar content of a cluster and what star formation efficiency is required 
for a bound object to remain. 

Using supercomputers N = 10 7 or more is possible, giving the resolution to model fully 
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the dynamical and chemical evolution of a dwarf galaxy. In addition feedback can be placed 
into simulations of the formation of the first cosmological objects. 
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1 2 3 4 5 6 

Time/Myr 

Fig. 1. — The change in the total energy with time for the ten SNRs detailed in table [l] 
normalised to the total energy within the simulation volume at 1 Myr in units of 10 51 ergs. 
The 'lost' energy has been radiated away. 
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1 r 




Radius/pc 

Fig. 2. — The energy balance for the SNRs in runs 1002 (cross), 1017 (square) and 1019 
(diamond) at 1 Myr. The solid lines show the sum of the thermal and kinetic energy, the 
dashed lines show the thermal and the dash-dot line the kinetic energy within each radius. 
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Fig. 3. — The evolution of the shell radius with time for the ten convergence runs detailed in 
table [I]. The fits to different evolutionarystages are explained in the text. The rough region 
of transition fromthe radiative pressure driven snowplough phase to the constant expansion 
phase is indicated by the dashed line at 1.5 Myr. 
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Fig. 4. — The evolution of the shell radius for 8 different ambient density ISMs all at 10 4 
K detailed in table fj. The hydrogen atom density per cm 3 is marked by each track. The 
approximate time and radius at which the shell velocity falls to the ISM sound speed (and 
the tracks become curved) is marked by the dashed line. 
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Fig. 5. — The total (solid line), thermal (dashed line) and kinetic (dot-dash line) energies 
within each radius at a shell radius of ~ 2 x tpds the radial distribution being given in terms 
of r/rpDs for three different densities of no = 0.01 cm~ 3 (triangles), no = 0.5 cm -3 (stars) 
and n = 100 cm~ 3 (circles). 
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Fig. 6. — The radial energy profiles of runs 1003 (stars) and 1018 (triangles) at 1 Myr 
compared to that of a run started at r = 2 x tpds (circles). The shell radius at this time 
is ~ 55 pc. The energies shown are the total (full line), thermal (dashed line) and kinetic 
(dot-dash line) within each radius. 
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Fig. 7. — The evolution of the shock velocity with radius for the standard M 5 = 1, R = 3.5 
cloud. 
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Fig. 8. — The mass loss rate across the sides of the simulation cube with time for a 
M = 1O 5 M , R = 3.5 pc, r max = 20 pc Plummer cloud with a central supernovae. 
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Fig. 9. — The minimum shock velocity v min in km s 1 against the Plummer cloud 
characteristics Rpf 2 /M^ 77 from eqn. 15. 
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Fig. 10. — The fractional mass lost against the ratio of the minimum velocity to escape 
velocity for the Plummer clouds listed in table [| A fit from eqn. is added for clouds with 
0.2u esc < f min < t> esc (clusters with partial mass loss - see text). 



-42- 



0.03 - 



CD 
0) 

0.02 

if) 

o 



0.01 




0.06 



0.14 



T min / 10 e rgs 



Fig. 11. — The feedback energy E\ osXl against the minimum kinetic energy parameter 
Tmm = O.OOIM5 (t> min /km s -1 ). Linear fits of the same slope are fitted from eqns. |T7| & 
TBI are marked by dashed lines. 
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Fig. 12. — The final-to-initial potential energy ratio against the ratio of the minimum velocity 
to escape velocity for the Plummer clouds listed in table [| The linear relationship marked 
by the dashed line is from eqn. ITS. 
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Fig. 13. — The density distribution after 20 Myr for the 1O 5 M Plummer cloud from fig. [8[ 
The solid line is a Plummer model with M = 6.36 x 1O 4 M and R = 7.8 pc 



-45- 



Table 1. Parameters of the ten uq = 0.5 cm -3 , ( m = 1 convergence test runs. The table 
gives the run identification number, box size (in pc), softening (in pc), number of particles 
N, the end time of the simulation (in Myr), the mass of an individual particle (in M & ) and 
the symbol used to represent each run in figs. and |j. 



RUN 


box 


softening 


N 


End time 


Particle mass 


Symbol 




pc 


pc 




Myr 


M Q 




1001 


300 


27 


32 3 


3 


11.1 


dot 


1002 


300 


23 


2 x 32 3 


3 


5.53 


cross 


1003 


300 


18 


4 x 32 3 


3 


2.77 


star 


1004 


300 


5 


64 3 


2 


1.38 


circle 


1005 


300 


5 


2 x 32 3 


3 


5.53 


'x' 


1017 


200 


11 


4 x 32 3 


2 


0.82 


square 


1018 


400 


20 


4 x 32 3 


5 


6.56 


triangle 


1019 


500 


28 


4 x 32 3 


6 


12.8 


diamond 


1020 


500 


22 


64 3 


6 


6.40 


open star 


1021 


200 


9 


64 3 


2 


0.41 


open cross 
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Table 2. The density of the surrounding ISM, the initial radius of the PDS and the 
approximate time and radius at which the shell velocity falls below the ambient medium's 
sound speed. The evolution of the shell radiuswith time for the different ambient densities 

used here is plotted in fig. f|. 



RUN 


n 


?~PDS 


^sound 


^*sound 




cm -3 


pc 


Myr 


pc 


1006 


0.01 


101 


7.2 


360 


1007 


0.1 


38 


3.0 


150 


1008 


0.2 


28 


2.2 


120 


1003 


0.5 


18 


1.6 


84 


1009 


1 


14 


1.3 


66 


1010 


5 


7.0 


0.70 


37 


1011 


10 


5.2 


0.45 


26 


1012 


100 


1.9 


0.21 


11 
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Table 3. The cloud models and feedback parameters of the runs in Section 3. The 
Plummer parameters mass M 5 and radius R pc for each cloud are followed by the minumum 
velocity t> min of the shell in that cloud, the fractional mass loss M\ ost /M, the energy 
feedback E\ ost and the loss of potential energy f2i ost from that cloud. 



M 


R 


^min 


M lost /M 


^lost 


^lost 


1O 5 M 


pc 


km S _1 




10 48 ergs 


10 48 ergs 


0.5 


2 


14.7 


1.00 


19 


32 


0.5 


3.5 


16.2 


1.00 


35 


18 


0.75 


2 


9.9 


0.17 


7.9 


45 


0.75 


3.5 


12.3 


1.00 


26 


40 


0.9 


2.5 


9.6 


0.16 


8.2 


50 




2 


8.3 


0.076 


5.1 


50 




2.25 


7.8 


0.080 


4.7 


45 




2.25 


8.2 


0.092 


5.7 


49 




2.5 


9.3 


0.13 


7.4 


53 




3.5 


10.2 


0.33 


15 


57 




3.5 


10.0 


0.32 


15 


56 




3.5 


10.2 


0.33 


15 


57 




5 


10.9 


1.00 


28 


49 


1.2 


3.5 


8.4 


0.16 


9.5 


60 


2.5 


3.5 


4.9 


0.014 


1.6 


56 


5 


3.5 


2.4 








34 


5 


5 


3.6 


0.004 


0.5 


62 



